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We report on a numerical study of real-time dynamics of electromagnetically interacting chirally 
imbalanced lattice Dirac fermions within the classical statistical field theory approach. Namely, we 
perform exact simulations of the real-time quantum evolution of fermionic fields coupled to classical 
electromagnetic fields, which are in turn coupled to the vacuum expectation value of the fermionic 
electric current. We use Wilson-Dirac Hamiltonian for fermions, and non-compact action for the 
gauge field. In general, we observe that the backreaction of fermions on the electromagnetic field 
prevents the system from acquiring chirality imbalance. In the case of chirality pumping in parallel 
electric and magnetic fields, electric field is screened by the produced on-shell fermions and the 
accumulation of chirality is hence stopped. In the case of evolution with initially present chirality 
imbalance, axial charge tends to transform to helicity of electromagnetic field. By performing 
simulations on large lattices we show that in most cases this decay process is accompanied by 
the inverse cascade phenomenon which transfers energy from short-wavelength to long-wavelength 
electromagnetic fields. In some simulations, however, we observe a very clear signature of inverse 
cascade for the helical magnetic fields which is not accompanied by the axial charge decay. This 
suggests that the relation between inverse cascade and axial charge decay is not as straightforward 
as predicted by the simplest form of anomalous Maxwell equations. 

PACS numbers: 12.38.Aw, ll.15.Tk 


I. INTRODUCTION AND BRIEF SUMMARY 

Over the past few decades, real-time instability of the 
system of chiral fermions coupled to dynamical gauge 
fields has been attracting a lot of attention in various 
fields of physics, ranging from astrophysics to condensed 
matter physics. This instability manifests itself in the de¬ 
cay of the initial imbalance between the densities of the 
left- and right-handed fermions at the expense of the gen¬ 
eration of magnetic fields with nonzero magnetic helicity 
(or, in other words, winding number of magnetic flux 
lines). In the astrophysical context, the phenomenon of 
chiral plasma instability is actively discussed as the mech¬ 
anism responsible for the generation and enhancement of 
primordial magnetic fields [1-5] as well as for the transfer 
of magnetic field energy from short to cosmological scales 
[3, 6, 7]. 

In the context of condensed matter physics, chiral 
plasma instability was initially discussed and experimen¬ 
tally detected as the so-called helical instability of liquid 
^He [8]. It was also considered as a mechanism of spon¬ 
taneous magnetization of topological magnetic insulators 
[9]. In experiments in which chirally imbalanced Weyl 
semimetal states are created from Dirac semimetals by 
applying parallel electric and magnetic fields [10-12] chi- 
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ral plasma instability might manifest itself in the spon¬ 
taneous emission of circularly polarized terahertz-range 
electromagnetic radiation [13]. 

In heavy-ion collisions, chiral plasma instability might 
lead to enhanced emission of circularly polarized soft pho¬ 
tons [13]. It should be also important for the correct esti¬ 
mate of the lifetimes of chirality imbalance and magnetic 
fields [14, 15]. However, the estimates of [16] suggest that 
in heavy-ion collisions the volume and the lifetime of the 
quark-gluon plasma might be too small for the instability 
to develop. 

The origin of this instability of chirally imbalanced 
Dirac fermions is the Chiral Magnetic Effect (CME) 
[17, 18] - electric current flowing parallel to the magnetic 
field in the presence of chirality imbalance. Within the 
linear response approximation the contribution of CME 
to the electric current is 

jCME = CrcME B. (1) 

The commonly quoted value for the chiral magnetic con¬ 
ductivity acME is <JCME = 2 ^) where / j,a is the so-called 
chiral chemical potential which parameterizes the differ¬ 
ence between the Fermi levels of right- and left-handed 
fermions and hence also the total axial charge of the sys¬ 
tem. The value of acME, however, strongly depends on 
frequency w and wave vector k of electromagnetic field, 
and, in the limit of constant and homogeneous magnetic 
field, on the way in which the limits w —t 0 and fc —>■ 0 
are taken [18-24]. 

In order to see how the CME current (1) can lead to in- 
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stability, one can insert it into the classical Maxwell equa¬ 
tions, along with the conventional Ohmic current j = aE, 
where a is the electric conductivity. Assuming the unbro¬ 
ken translational invariance both in time and space, we 
can write these so-called anomalous Maxwell equations 
[1, 3, 4, 13, 14, 16, 25-27] in frequency-momentum space 
as 

iwB = —ik X E, iwE = ik x B — a E — acME B. (2) 

From these equations we find the following four-branch 
dispersion relation for transversely polarized plane waves 
with the wave vector k = k [16]: 

icT / „ cr^ , . 

= y + sy -k rcrcMfifc - (3) 

where s = ±1 and r = ±1 label different branches of 
the dispersion relation. The corresponding polarization 
vectors = 2“^/^ (1, —zr, 0) for the electric field E cor¬ 
respond to circularly polarized waves with opposite he- 
licities (handedness) for opposite r. 

While for nonzero electric conductivity a the imagi¬ 
nary part of w in (3) is always positive and hence corre¬ 
sponds to decaying plane waves, nonzero chiral magnetic 
conductivity can also lead to exponentially growing solu¬ 
tions if the absolute value of the wave vector k is smaller 
than acME- From (3) it is also easy to see that for a 
given wave vector k, only one of two helical mode will 
exhibit exponential growth. For example, for fiA > 0 
(and hence Qa > 0 and ucme > 0 ) and <tcme > k > 0 
the exponentially growing solution has the form 

El = /e''* cos (kxs ), E 2 = —/e'^* sin {kxs ), 

Bi = —/ —e”* cos (kxs ), B 2 = /—sin (kx ^), 

K K 

A 3 = 0, A 3 = 0, (4) 

where / is an arbitrary amplitude and k = —iw = 
— f -I- ^ — k"^ + acMEk. It is important to stress that 

since this solution grows monotonously in time, here we 
use the terms “circular polarization”, “handedness” and 
“helicity” to describe the rotation of the vectors E and 
B along the X 3 axis, rather than in time. The growth 
of long-wavelength electromagnetic waves and the de¬ 
cay of short-wavelength waves predicted by the anoma¬ 
lous Maxwell equations (2) is a novel mechanism for 
the inverse cascade in relativistic magnetohydrodynam¬ 
ics [3, 6 ], which transfers energy from long- to short- 
wavelength helical magnetic fields. 

The fact that the exponentially growing solution (4) 
has the helical structure of the form (4) also suggests the 
mechanism which can stop the growth of electromagnetic 
field at later times. Namely, let us recall that for massless 
chiral fermions the time evolution of the axial charge is 
governed by the anomaly equation: 

OiQa = ^ J d^xE ■ B, (5) 


where the axial charge Qa = Qr — Ql is defined as 
the difference between the charges Qr and Ql of the 
right- and left-handed fermions, g is the electromagnetic 
coupling constant and we have integrated over space to 
get rid of the spatial divergence of the axial current. For 
simplicity, in this paper we consider only a single flavor 
of Dirac fermions with electromagnetic coupling g = 1. 

For the exponentially growing solution (4) the product 
E ■ B is negative: E ■ B = _j 2 ^g 2 Kt anomaly 

equation (5) then dictates that the time derivative OiQa 
of the axial charge is negative. Since we have assumed 
Qa > 0 , /za > 0 , we see that the growing helical solu¬ 
tion (4) will result in the decrease of Qa and hence of 
fiA- This depletion of chirality imbalance should even¬ 
tually suppress the chiral magnetic conductivity in ( 1 ) 
and hence slow down or stop completely the exponential 
growth in (4). 

However, the above analysis of the chiral plasma insta¬ 
bility, which follows [3, 4, 13, 14, 16, 25-28], essentially 
relies on an assumption that the electric current takes 
the form j = aE -|- acMEB with constant ohmic and 
chiral magnetic conductivities. In reality, both cr and 
crcME depend on the frequency and wave vector of elec¬ 
tromagnetic field in a nontrivial way [18-24]. One can 
also expect a strong dependence of a and cjcme on the 
spatial and temporal modulation of the axial charge den¬ 
sity, which will in general appear at late evolution times 
[27]. Moreover, as the instability might lead to quite 
large strengths of electric and magnetic fields, nonlinear 
effects beyond the linear response result ( 1 ) might be¬ 
come important. Using linear response approximation to 
describe the interactions between the fermions and the 
electromagnetic fields is in fact similar to the Lyapunov 
analysis of the full quantum evolution, which is in general 
nonlinear. What concerns inter-fermion interactions, so 
far they were taken into account only indirectly, by using 
the relaxation time approximation [15] or the decoher¬ 
ence of the fermionic wave functions [29]. A consistent 
inclusion of all these effects in the anomalous Maxwell 
equations ( 2 ) would be certainly difficult with approaches 
based e.g. on the chiral kinetic theory [25-27, 30], chi¬ 
ral hydrodynamics [ 6 , 31] or the Langevin-type effective 
theory [32]. 

This situation clearly calls for a more first-principle de¬ 
scription of the real-time dynamics of chirally imbalanced 
plasma which would overcome these limitations and ap¬ 
proximations. In this paper, we report on the numerical 
study of the real-time chiral plasma instability within the 
framework of the so-called classical statistical field theory 
(CSFT) [33-38], which captures the first nontrivial order 
of the expansion of the full quantum evolution operator 
in powers of the Planck constant. CSFT is currently the 
state-of-the art method for numerical simulations of real¬ 
time quantum evolution. The CSFT approximation is 
justifiable as long as the characteristic occupation num¬ 
bers of the physically relevant gauge field modes are large. 
That is, the dynamics of the gauge fields should be almost 
classical. On the other hand, the real-time dynamics of 
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FIG. 1: Two ways of introducing initial chiral imbalance for 
the many-body Dirac Hamiltonian. On the left: by intro¬ 
ducing the chiral chemical potential /ia in the single-particle 
Dirac Hamiltonian. On the right: by filling more right-handed 
eigenstates and less left-handed eigenstates (or vice versa). 


fermions is exact in CSFT. Taking into account that in 
all previous studies gauge fields were also treated classi¬ 
cally, the applicability of the CSFT approach is obviously 
wider than that of the previously used approaches. An 
approach very similar to CSFT has been recently used 
in [39] to study the real-time dynamics of CME. How¬ 
ever, in this work the back-reaction of fermions on the 
electromagnetic field, which is the origin of the chiral 
plasma instability, was not taken into account. The real¬ 
time dynamics of axial charge was also studied in 1 -b 1- 
dimensional Abelian Higgs model in the pioneering work 
[33]. 

Our studies are based on the non-compact formulation 
of lattice quantum electrodynamics, which avoids poten¬ 
tial problems with monopole condensation in the strong¬ 
coupling phase [40]. For fermions, we use the massless 
Wilson-Dirac Hamiltonian which has a low-energy chiral 
symmetry. At sufficiently high energies, this symmetry is 
broken due to the Wilson term. In the condensed matter 
context, this breaking is a natural feature of any model 
description of Dirac and Weyl semimetals [41-43]. In 
Section HI we demonstrate that the effect of this explicit 
breaking is, however, not very large (see also [44]). There¬ 
fore we hope that our results should be also at least qual¬ 
itatively relevant in the context of high-energy physics, 
where the chiral symmetry is exact at the level of the 
Lagrangian, or tends to be exact at sufficiently high en¬ 
ergies. 

In order to introduce the initial chirality imbalance, we 
have started the simulations with a state in which more 
right-handed eigenstates and less left-handed eigenstates 
are filled, as depicted on the right panel of Fig. 1. Such a 
state is an excited state of the many-body Dirac Hamil¬ 
tonian, even in the absence of electromagnetic fields. 
It is an idealized description of the result of “chirality 
pumping” process in parallel electric and magnetic fields 
[10, 45] or in intense circularly polarized laser beams. 

Alternatively, we have also considered the introduction 
of the chiral chemical potential into the single-particle 
Dirac Hamiltonian, which changes the energies of the 
right- and the left-handed Dirac points (see left panel 


of Fig. 1). Such initial state has nonzero axial charge but 
is still the ground state of the many-body Hamiltonian 
in the absence of interactions with electromagnetic fields. 
In simulations which started from this state we have not 
found any signatures of instability or the transfer of he- 
licity from fermions to electromagnetic fields. The ax¬ 
ial charge density exhibited only small fluctuations on 
top of the large mean value. Presumably, the reason for 
such a behavior is that nonzero chiral chemical potential 
corresponds to the physical situation in which our sys¬ 
tem is connected to an infinite reservoir of axial charge, 
which is capable of maintaining its initial value at a con¬ 
stant level. Since the anomaly equation (5) holds also at 
nonzero chiral chemical potential, this also implies that 
the magnetic helicity can only exhibit small fluctuations, 
possibly related to the violation of the anomaly equation 
due to lattice artifacts. Since such behavior is not really 
interesting, we do not discuss this setup in what follows. 

The structure of this paper is the following: in Sec¬ 
tion H we start with a brief summary of the details of 
our numerical CSFT algorithm. In Appendix A we pro¬ 
vide a more detailed derivation of this algorithm with a 
bias towards non-relativistic field theories and condensed 
matter systems, which might hopefully complement the 
existing literature on CSFT (see e.g. [33, 38] for deriva¬ 
tions which are more in the spirit of relativistic quan¬ 
tum field theory). In this Appendix we also demonstrate 
explicitly the absence of any nontrivial Jacobian in the 
integration measure in the CSFT algorithm, and discuss 
some practical aspects of our CSFT simulations on paral¬ 
lel computers. In Section HI we present the results of the 
simulations of the chirality pumping process. First, we 
consider chirality pumping in external parallel electric 
and magnetic fields in the absence of backreaction and 
verify the validity of the anomaly equation (5) in our nu¬ 
merical setup. After that we consider the effect of back- 
reaction of fermionic current on the chirality pumping 
process and demonstrate that the dynamical screening 
of the external electric field prevents the system from ac¬ 
quiring large axial charge density at late evolution times. 

In Section IV, we consider the decay of the initial chi¬ 
ral imbalance and the generation of electromagnetic fields 
with nonzero helicity. In order to trigger the decay, we 
start simulations with several initially excited modes of 
electromagnetic field. Following the energies of helical 
modes in momentum space, we demonstrate that only 
long-wavelength modes of definite helicity grow and all 
other modes decay with time. This is a direct numerical 
evidence of the inverse cascade phenomenon [3, 6, 13, 46] 
due to chiral plasma instability. We find, however, that 
the dependence of the strength of the inverse cascade 
on the initial conditions and parameters of the simu¬ 
lations is significantly more complex than predicted by 
the anomalous Maxwell equations (2). In particular, in 
the simulations which exhibit most rapid growth of he¬ 
lical magnetic fields the axial charge does not decay at 
all. Correspondingly, the mechanism which stops the in¬ 
verse cascade in our simulations is not related to the axial 
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charge decay, again in contrast to the expectations based 
on the equations (2) [3, 4, 6 , 13-15]. Rather, we observe 
that in most simulations which do exhibit axial charge 
decay the inverse cascade emerges for the electric rather 
than for magnetic helds. In Section V we conclude with 
a general discussion of our results and an outlook. 


II. CLASSICAL STATISTICAL FIELD THEORY 
APPROXIMATION TO REAL-TIME 
EVOLUTION 

We consider the many-body fermionic Hamiltonian 
coupled to dynamical non-compact electromagnetic helds 
on the lattice, so that the full Hamiltonian H of our sys¬ 
tem is H = Hp + HpM- The fermionic Hamiltonian Hp 
reads 

Hp = ( 6 ) 

x,y 


The operators of the electric held Ex,i and the vector 
potential are canonically conjugate and satisfy the 


commutation relations 


Exji ; AyJ 


= —iSxySij. We im¬ 


pose periodic boundary conditions in all spatial direc¬ 
tions both for the gauge and the fermionic helds. 

In our CSFT algorithm, described in detail in Ap¬ 
pendix A, we numerically solve the classical equations of 
motion of the electromagnetic held with the Hamiltonian 
(9) 


d^Ax^t {t) = -Jx,i (t) - {jx,z {t )) - 

~ {Px,ij {t) — Fx-epij {t)) , ( 11 ) 

3 


where the initial value of the time derivative 
dtAx,i{t)\^^Q is the initial value of the electric held 
Ex,i ( 0 ) and {jx,i (t) ) is the vacuum expectation value of 
the fermionic electric current, which can be calculated 
as 


where the labels a;, y denote the sites of the three- 
dimensional cubic lattice, ipx are the spinor-valued 
fermionic creation and annihilation operators which sat¬ 
isfy the anti-commutation relation = 5xy and 

hx^y is the massless single-particle Wilson-Dirac Hamil¬ 
tonian with the Wilson coefficient r = 1: 

3 

^x,y — 3vp j35x^y ^ ^ 4- O-i) e ^ 5y^x+ei T 

i=l 

i=l 

Here Ax^i is the vector potential of the lattice gauge held, 
ei denotes the unit lattice vector in the direction i, vp 
is the Fermi velocity, jS and are the Dirac /3- and a- 
matrices and 75 is the generator of chiral rotations: 



where Ui are the Pauli matrices. In (7) we have assumed 
that the lattice spacing is unity. Thus in what follows 
all dimensionful quantities are expressed in units of the 
lattice spacing. 

The lattice Hamiltonian Hem of electromagnetic held 
is the straightforward lattice discretization of the corre¬ 
sponding continuum Hamiltonian: 


Hem = E E 

X 2 = 1 


(A ■ ^ ■ ■ 

^ + 1 :^ 

\ 



, ( 9 ) 


(3 x,i {t)) = Tr (po u (0, t) jx,i (0, t)) , (12) 

where jx,i = is the single-particle operator of elec¬ 
tric current, u(0,t) is the quantum evolution operator 
dehned by the single-particle Schrbdinger equation 

dtu(0,t) = ihlAx,i(t)]u(0,t), u(0,0) = I, (13) 

and po is the initial density matrix which characterizes 
the initial occupation numbers Ua of single-particle states 

lA)'- 

Po = E 

a 

In our case, \ipa) are the eigenstates of the single-particle 
Hamiltonian (7). If some occupation numbers are exactly 
zero (which can be the case at zero temperature), some 
components of the quantum evolution operator u com¬ 
pletely decouple and can be discarded in the solution of 
the equation (13). This can be used to speed up the 
algorithm, typically by a factor of two (for a standard 
zero-temperature Fermi distribution). The expectation 
value {jx^i (t)) in ( 11 ) describes the effect of backreac- 
tion of fermions on the electromagnetic fields. 

We thus have a closed set of equations (11), (12) and 
(13), which allows to evolve the fermionic quantum states 
and the classical electromagnetic fields in a self-consistent 
way. One can also explicitly check that this evolution 
conserves the total energy Hem + {Hp) of electromag¬ 
netic field and fermions up to the work done by the ex¬ 
ternal current (t): 


where Jx,i (t) is the external current (which is required, 
e.g., to switch the external electric and magnetic fields on 
and off) and the operator of the magnetic field strength 
tensor Fx,ij is defined in terms of the finite differences of 
the vector potential operator Ax,i as 

^x,ij — Ax^i -f Ax-\-ei,j Ax-\-ej ,i Axg- (f9) 


dt [Hem + (ilF)) = - E (^) (^) ^ 


Hem = 2 E ( AtAx,i)^ + E ^x,ij 

x,i \ j 

{Hp) =Tt {pou{0,t) huA0,t)) ■ ( 15 ) 
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We have solved the evolution equations (11) and (13) 
using the leap-frog integrator, which slightly violates the 
conservation of energy (15). At sufficiently small time 
step this violation is completely under numerical control, 
see Fig. 14 in the Appendix A. 

In the CSFT approach, one can also partially take 
into account the quantum fluctuations of the electromag¬ 
netic fields, encoded in the nontrivial Wigner transform 
Pem i^o, Eg) of the initial density matrix pem, where 
Aq = A{t = 0) and Eq = E {t = 0) are the initial values 
of electric and magnetic fields. To this end one should 
additionally average all observables over Aq and Eq, sam¬ 
pled with the probability pem (Aq, Ao)[58]. However, in 
our work we have not taken these initial quantum fluc¬ 
tuations into account for the following reasons. First, in 
the case of anomaly equation (5) with massless fermions 
we have found that the effect of initial fluctuations of 
electromagnetic fields is much more significant than in 
the case of e.g. Schwinger pair production [38] [59], and 
hence much more samples of the initial fields are required 
to reach acceptable statistical errors. Partially this can 
be explained by the large value of the electromagnetic 
coupling constant, which was g = 1.0 in most of our sim¬ 
ulations. We expect that the role of initial fluctuations 
will be smaller for a smaller value of g, say, 5 = 0.1. In 
the latter case, however, the characteristic time scale of 
the chiral plasma instability increases significantly above 
our current simulation times. 

In addition, taking into account initial quantum fluctu¬ 
ations of electromagnetic fields makes it impossible to as¬ 
sume spatial homogeneity of electromagnetic fields along 
some of the lattice directions, which is essential to speed 
up the CSFT simulations at large lattice sizes. 

Thus while the effect of quantum fluctuations on the 
chiral plasma instability might be potentially very sig¬ 
nificant and interesting, we cannot study it with our 
presently available computational resources and leave it 
for the future work. In this work, we avoid the statis¬ 
tical averaging over the initial values of the fields Eq, 
Aq by using the very simple form of the initial density 
matrix pem [Aq, Eq) which is just a delta-function on 
some particular, specifically chosen initial values. We 
thus completely neglect the quantum fluctuations of the 
electromagnetic fields. Nevertheless, this approximation 
is still certainly wider than the chiral kinetic theory or 
hydrodynamical approximation. 


III. CHIRALITY PUMPING IN PARALLEL 
ELECTRIC AND MAGNETIC FIELDS 

In this Section we study the real-time evolution of the 
axial charge Qa in the background of constant parallel 
external electric and magnetic fields. In the absence of 
backreaction, such setup provides a direct check of how 
well the anomaly equation (5) holds for the Wilson-Dirac 
Hamiltonian with inexact chiral symmetry [44] , which we 
further use to study the chiral plasma instability in Sec¬ 


tion IV. The effect of backreaction is also interesting since 
the anomaly equation (5) is known to receive nontriv¬ 
ial corrections if the electromagnetic fields are dynamical 
[47-49]. 

In order to induce the constant external electric field 
E = E € 3 , we switch on the external current of the form 
J^x,i (t) = Si,3 E t. Constant external magnetic field is 
induced by the static circular external current flowing 
around the plaquettes with xi = Li — 1 , X2 = L2 — 1 
for all cca = 0 ... L 3 — 1 , where Li, L 2 and L 3 are lattice 
sizes. This static current is like a thin solenoid piercing a 
stack of lattice plaquettes, with the field strength being 
equal to B outside of solenoid and B — BL 1 L 2 inside it. 
Lattice fermions, however, acquire only the Aharonov- 
Bohm phase e*®® when encircling such a solenoid if one 
imposes the flux quantization condition 

gBLiL 2 = 27r$, $ S Z. (16) 

The total external current which we insert in the equa¬ 
tions (If) is the sum of the two currents which create 
constant electric and magnetic fields. 

For the initial state of the fermionic fields, we use 
the eigenstates of the Wilson-Dirac Hamiltonian h [Ag], 
where Aq is the initial gauge field configuration with con¬ 
stant magnetic field B as described above. In this work 
we consider only the limit of zero temperature, corre¬ 
spondingly, only the states with negative energies are 
initially occupied. 

The Wilson-Dirac Hamiltonian which we use in our 
simulations does not have exact chiral symmetry, and 
there is no uniquely defined axial charge operator which 
would exactly satisfy the anomaly equation (5) and com¬ 
mute with the Hamiltonian. Rather, the anomaly equa¬ 
tion (5) can only hold approximately, in the limit of large 
lattice volume and sufficiently smooth, slowly changing 
and small gauge fields [50, 51]. We thus use the simplest 
possible definitions of the operators of the axial charge 
density qax and the total axial charge Qa- 

qAx=i’\lb^x, QA = ^qAx- (17) 

X 

The time-dependent expectation value of the axial charge 
density is calculated similarly to the expectation value of 
electric current in ( 12 ): 

{qAx{t)) = Tr (pou(0,t) (0,t)) , (18) 

where Px is the single-particle projector on a single lattice 
site x: [Px]xiX 2 ~ ^xix^xx 2 - 

First we neglect the backreaction of fermionic electric 
current on the electromagnetic field and measure the time 
dependence of the axial charge in constant parallel exter¬ 
nal electric and magnetic fields. The results are shown 
on the left panel of Fig. 2 for the 10 x 10 x 32 lattice with 
$ = 1 quantum of magnetic field flux. One can see that 
Qa grows linearly with time until it reaches some maxi¬ 
mal value QaIV ~ 0.006, where V = L 1 L 2 L 3 is the total 
number of lattice sites (lattice volume). After that, Qa 
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FIG. 2: Chirality pumping without backreaction for various external electric helds and external magnetic field with fiux $ = 1 
on the 10 X 10 X 32 lattice. On the left: time dependence of the axial charge Qa and its linear fits at early times. On the 
right: dependence of the slope of these fits on the external electric field and the linear fit of this dependence. Here and in what 
follows time is expressed in units of lattice spacing. 


decreases again. This decrease is a lattice artifact related 
to the fact that the characteristic momentum p ^ Et oi 
fermions accelerated by an electric field E approaches 
the UV cutoff set by the compact size of the lattice mo¬ 
mentum space G [—tt ... tt]. Due to the periodicity of 
lattice momentum space, at large time scales the behav¬ 
ior of the axial charge (in the absence of backreaction) is 
well described by Qa ~ sin (M/2). There are also some 
short-time fluctuations on top of the clearly visible linear 
growth at early times. 

In order to estimate the linear growth rate at early 
times, we perform the linear fit of the form Qa (t) jV = 
a (E) t in the range t G [0... 50] for £1 = 0.01 and 
E = 0.02 and in the range t G [0... 30] for other val¬ 
ues of E. The dependence of the coefficient a {E) on 
the electric field is shown on the right panel of Fig. 2. 
Again, this dependence is linear with a good precision, 
and we perform another linear fit a (E) = CE, where C 
corresponds to the anomaly coefficient relating OiQa and 
/ (PxE-B in (5). On Fig. 3 we show the dependence of C 
on the size of the lattice (in the directions perpendicular 
to the magnetic field). One can see how C approaches 
the value C = in the limit of large lattices, in agree¬ 
ment with the anomaly equation (5). Let us also note 
that for larger number of flux quanta one can perform 
a similar fitting procedure. However, finite-volume ar¬ 
tifacts in C are significantly larger for larger magnetic 
fluxes. For this reason, in this work we have only used 
external magnetic field with one flux quantum. 

It is also interesting to check how the axial charge de¬ 
pends on time after the external electric field is switched 
off and the external magnetic field remains constant 
(which we believe to be a more realistic experimental 
setup than the simultaneous switching off of all fields). 
The time dependence of the axial charge for such a sit¬ 
uation is shown on the left plot of Fig. 4 (green line). 
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FIG. 3: Dependence of the lattice anomaly coefficient 
{BiQa = C J (ExE ■ B) on the transverse lattice size. Lat¬ 
tice size in the direction parallel to the magnetic field is fixed 
to La = 32. 


External electric field is switched off at the time t = 50. 
One can see that starting from this moment of time the 
axial charge exhibits only some small-scale fluctuations 
around the nonzero mean value. This demonstrates that 
the effect of explicit chiral symmetry breaking due to the 
Wilson term in the Hamiltonian (7) is rather small for 
such simulation parameters, and the total axial charge is 
almost a conserved quantity. 

After establishing the validity of the anomaly equa¬ 
tion (5) in our simulation setup, we study the effect of 
backreaction of dynamical electromagnetic fields on the 
chirality pumping process. Technically, the backreaction 
is taken into account by inserting the expectation value 
of the fermionic electric current {jx,i) into the Maxwell 
equations for the electromagnetic field. We now consider 
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FIG. 4: A comparison of chirality pumping processes with and without backreaction of the fermionic electric current on the 
electromagnetic field. On the left: time dependence of the axial charge. On the right: time dependence of the component of 
the volume-averaged electric held parallel to the magnetic held. Lattice size is 10 x 10 x 32, the hux of external magnetic held 
is <E> = 1, external electric held is 0.01. 


the situation in which the external electric and magnetic 
fields are switched on permanently. As we will see, in 
simulations with backreaction switching off the electric 
field at sufficiently late times anyway does not affect the 
evolution significantly due to screening by the dynam¬ 
ically generated electric field. Time dependence of the 
axial charge Qa (i) for simulation with backreaction is 
shown on Fig. 4. For comparison, on the same Figure we 
also show Qa (t) for simulations without backreaction, 
where the electric field is permanent or switched off at 
t = 50. 

One can see that while at f < 30 Qa (i) grows approx¬ 
imately linearly with t both with and without backreac¬ 
tion, at later times backreaction leads to a rapid decay of 
Qa with subsequent fluctuations around zero. In order 
to understand the origin of this effect, remember that ax¬ 
ial anomaly can be also regarded as the Schwinger pair 
production in the effective 1 -I- 1-dimensional theory of 
fermions on the lowest Landau level. It is thus natural to 
expect that particle-anti-particle pairs produced by the 
external electric field will tend to screen this field, just 
as in the case of Schwinger effect in (3 -I- 1) dimensions. 
To check this conjecture, on the right panel of Fig. 4 we 
plot the volume-averaged electric field projected on the 
direction of the magnetic field. One can see that indeed 
it quite quickly decreases from the initial value E = 0.01, 
reaching zero at around t « 30 - exactly at the time at 
which the growth of the axial charge stops (see left panel 
of the same Figure). After that the electric field exhibits 
some fluctuations around zero with the amplitude which 
is approximately five times smaller than the initial field 
value. We thus conclude that the effect of backreaction 
on the chirality pumping is to stop the growth of the ax¬ 
ial charge by screening the external electric field down to 
zero. 


IV. CHIRAL PLASMA INSTABILITY AND 
DECAY OF AXIAL CHARGE 

In this Section, we consider a situation in which some 
initial chiral imbalance is already created e.g. by chiral¬ 
ity pumping, and the parallel electric and magnetic fields 
are adiabatically switched off while keeping nonzero value 
of the total axial charge Qa and hence the chiral chem¬ 
ical potential fj-A- In this setup we would like to study 
the existence and the late-time evolution of the exponen¬ 
tially growing solutions (4) of the anomalous Maxwell 
equations ( 2 ), as well as the associated inverse cascade of 
energy of helical electromagnetic fields. 

In order to implement the initial chirality imbalance as 
discussed in the introductory Section I (see right panel 
of Fig. 1), we divide all the eigenstates of the single¬ 
particle Wilson-Dirac Hamiltonian /i[Ao], where Aq is 
the initial value of the vector potential, into the positive 
chirality states with {tpa \ 75 IV’a) > 0 and the negative chi¬ 
rality states with {-ipal 75 |^/’a) < 0. For positive chirality 
states we fill all the levels with €a < fJ- a, and for negative 
chirality states - all the levels with Cq < —fJ,A- While 
the eigenstates of the Wilson-Dirac Hamiltonian are not 
in general the eigenstates of the 75 operator, for eigen¬ 
states with sufficiently small momenta our definition is 
maximally close to the notion of distinct Fermi levels of 
left- and right-handed continuum massless fermions. In 
practice, this definition is unambiguous as long as all the 
energy levels are non-degenerate. In the case of de¬ 
generate energy levels (as e.g. in the case of the Wilson- 
Dirac Hamiltonian with zero gauge fields or in the back¬ 
ground of a single plane wave), one can additionally ro¬ 
tate the eigenstates within the degenerate subspaces in 
order to maximize the absolute values of matrix elements 
(V'aldS IV'a). 

With lattice discretizations of the Dirac Hamiltonian 




it is not possible to have very large values of the chiral 
chemical potential /xa, since the dispersion relation at 
the Fermi energy > 1 starts deviating from the Dirac 

cone e (jtj = VF\k\ due to lattice artifacts. At = 2, 
the Fermi energy touches the lowest van Hove singularity 
(saddle point) of the dispersion relation, and the excita¬ 
tions around the Fermi surface no longer correspond to 
Dirac fermions. On the other hand, the solution (4) of 
the anomalous Maxwell equations (2) with the conven¬ 
tional value fJcME = 2 ^ of chiral magnetic conduc¬ 
tivity suggests that the wave vectors at which the chiral 
plasma instability can occur are bounded by |A:| < 

On a finite spatial lattice of size L with periodic boundary 
conditions, the smallest nonzero value of |fc| is |fc| = 
which dictates the lower bound on the size of the lattice 
where the instability can be observed: 


L > 


47r^ 


(19) 


Thus it is advantageous to use large values of in or¬ 
der to reduce the lattice size used for simulations. Taking 
the moderate value ^a = 0.75, at which the dispersion 
relation is still linear with a good precision, we obtain 
L > 165. Performing simulations on an isotropic three- 
dimensional lattice of such a size would be a formidable 
numerical task. For this reason we have used the lat¬ 
tices with different sizes in different directions, so that 
the size L 3 in the direction 0:3 of electromagnetic wave 
propagation is much larger than the sizes Li = L 2 = 
in the transverse directions xi and X 2 - In addition, we 
have assumed that electromagnetic fields do not depend 
on the transverse coordinates xi and X 2 - This allows us 
to represent the single-particle evolution operator u ( 0 , t) 
in the block-diagonal form in the basis of plane waves 
propagating along xi and X 2 ^ which greatly reduces the 
dimensionality of the linear space on which the single¬ 
particle Schrodinger equation (13) should be solved. By 
comparing the results of simulations with — 20 and 
Lg = 40 at fixed L 3 = 200 (see Table I and Figs. 5, 
6 , 10 , 11 and 12 ) we have checked that the dependence 
on the transverse lattice size is rather weak. Let us also 
note that one of the reasons for not using the final state of 
chirality pumping process described in Section III for the 
study of chiral plasma instability is that in this case it is 
not possible to assume spatial homogeneity in transverse 
directions due to the breaking of translational invariance 
by external magnetic field [52]. 

While the fermionic initial state described above is an 
excited state which can spontaneously decay due to chiral 
plasma instability, in numerical simulations one always 
needs some small “seed” perturbation to start the decay 
process in a controllable way. For this reason we have 
started our simulations with a state in which also some fi¬ 
nite number n of electromagnetic field modes are excited. 
All of them are plane waves propagating along the lattice 
direction X 3 with the largest size L 3 , with a few small¬ 
est nonzero wave numbers km = , m = 1 .. .n and 


Set.No. 

L 3 

Lg 


n 

/ 

Vf 

Qa 4 

7ft 

'— 1 

1 

200 

20 

0.75 

10 

0.20 

1.00 

/ 

X 

/ 

2 

200 

40 

0.75 

10 

0.20 

1.00 

/ 

X 

/ 

3 

200 

20 

1.50 

10 

0.20 

1.00 

/ 

/ 

/ 

4 

200 

20 

0.75 

10 

0.05 

1.00 

X 

/ 

X 

5 

200 

20 

0.75 

4 

0.20 

1.00 

/ 

/ 

? 

6 

200 

20 

0.75 

4 

0.05 

1.00 

X 

/ 

? 

7 

200 

20 

1.50 

10 

0.05 

1.00 

X 

/ 

? 

8 

200 

20 

0.75 

10 

0.20 

0.75 

/ 

X 

/ 

9 

20 

20 

1.00 

1 

0.20 

1.00 

/ 

X 

X 


TABLE I: Summary of parameters and results of our simu¬ 
lations of chiral plasma instability. The column Qa 4 sum¬ 
marizes the decay of the axial charge and the columns t 
summarize the growth of the energies of long-wavelength elec¬ 
tric and magnetic fields. The symbols /, X and ? denote, 
respectively, the clearly visible growth, clearly visible absence 
of growth and intermediate situations for which it is difficult 
to make any conclusion within a finite simulation time. 


random linear polarizations. In order to facilitate the de¬ 
tection of the inverse cascade, we choose the amplitudes 
of all modes in such a way that their contributions to the 
total energy of electromagnetic field are equal. Thus the 
explicit form of our initial electromagnetic field configu¬ 
ration is 


A,., {t = 0)=Y, 


/ 


m—1 


’ (^m) 


i COS 0m) : 


Ea;,i (t = 0 ) = dt Aj,,, (t)lt^o = 
n 

= ^ / Sm{kmX3 + 4>m) , ( 20 ) 


where Umi are the random unit transverse polarization 
vectors which are chosen to coincide with one of the basis 
vectors ei, 6*2 with equal probability, G [ 0 , 27r] are the 

random phases and w (km) = y^4sin^ (%^) corresponds 
to the lattice dispersion relation for free massless fields 
on the lattice. 

In order to understand how the evolution process de¬ 
pends on various lattice parameters, we have performed 
simulations with 9 different parameter sets, which are 
summarized in Table I. We have varied both the trans¬ 
verse and the longitudinal lattice sizes, the initial elec¬ 
tromagnetic field amplitude /, the number n of initially 
excited electromagnetic field modes, the initial value of 
axial charge and the Fermi velocity vp. The parameter 
set No. 1 with L 3 = 200, Lg = 20, ^a = 0.75, n = 10, 
/ = 0.2 and = 1 is the “default” parameter set, and 
all other sets differ from it by a change in a few pa¬ 
rameters. Correspondingly, in what follows we label the 
data points on the plots which combine the results from 
several simulations by the number of parameter set (pre¬ 
ceded by the hash symbol #), in parentheses giving only 
those parameters which differ from the default ones. 
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FIG. 5: Time dependence of the total axial charge for different simulation parameters. On the right plot we compare simulations 
with different and linearly rescale qa (t) —>■ cqA {t) so that the initial values of cqA (0) agree for all simulations. In the plot 
labels, numbers preceded by the hash symbol # correspond to the numbers of parameter sets in Table I. In parentheses we give 
the values of only those parameters which differ from the default parameters (L 3 = 200, Ls = 20, fj,A ~ 0.75, n = 10, / = 0.2 
and vf = 1, parameter set No. 1). 


On Fig. 5 we show the time dependence of axial charge 
in simulations with parameters summarized in Table I. 
In simulations with initial amplitude of electromagnetic 
fields being equal to / = 0.2 the axial charge Qa de¬ 
cays with time. Interestingly, simulations with the small¬ 
est lattice size (parameter set No. 9) exhibit the fastest 
decay of Qa- On the other hand, with the initial am¬ 
plitude / = 0.05 the axial charge density exhibits only 
a rather small decrease at intermediate evolution times, 
subsequently followed by a slight increase. This non¬ 
trivial dependence on the electromagnetic held strength 
suggests that the dynamics of the decay process is more 
complicated than suggested by the anomalous Maxwell 
equations (2). It is interesting that the evolution of the 
axial charge seems to depend only weakly on simulation 
parameters other than the initial amplitude / and the 
longitudinal lattice size L 3 (through the value of the low¬ 
est wave number ). Even the dependence on the chiral 
chemical potential fj-A appears to be rather weak (after 
a trivial rescaling with respect to the initial value), see 
right plot on Fig. 5. The characteristic time scale for 
the evolution of the axial charge appear to be essentially 
larger than in the chirality pumping simulations in the 
previous Section. This difference can be qualitatively 
explained by much weaker field strengths in the simula¬ 
tions described in this Section. We also note that the 
initial values of the axial charge are roughly consistent 
with the continuum formula QaIV = 

V is the lattice volume. Deviations from this value can 
be explained, first, by the inclusion of the initial vector 
potential Aq in the initial Hamiltonian, and second, by 
the smaller value of chirality \{'tpa \ 75 \'>Pa)\ < 1 for high- 
energy eigenstates I'l/'a) of the Wilson-Dirac Hamiltonian 

( 7 ). 

A scaling analysis of the anomalous Maxwell equations 
suggests that at late evolution times the time dependence 



FIG. 6 : Universal late-time scaling Qa (t) ~ l/v^ of axial 
charge density in simulations with / = 0.2. Time depen¬ 
dence of (Qa ( 0 ) /Qa (t))^ in the second half of the total evo¬ 
lution time is htted, where appropriate, by linear functions 
(Qa (0) /Qa (t))^ = A + Bt. The hts are shown with dashed 
black lines. 


of the axial charge density approaches the simple power 
law [13, 46] 

Qa (t) ~ 1 /Vi. ( 21 ) 

In order to check this scaling, on Fig. 6 we plot the 
time dependence of the inverse square of the axial charge, 
which should approach the linear function I/ Q^ (t) ~ t 
according to (21). This asymptotic behavior indeed 
seems to emerge at late evolution times for simulations 
with Lg = 20, n = 10, / = 0.2 and ^a = 0.75, both with 
Vf = 1 and vp = 0.75 (parameter sets No. I and 8 ). The 
linear fits of (Qa ( 0 ) /Qa (t)) for these simulations are 
shown on Fig. 6 with dashed black lines. 
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We now check whether the decay of the axial charge is 
accompanied by the growth of the long-wavelength modes 
of electromagnetic field, as predicted by the anomalous 
Maxwell equations (2). To this end we perform the 
Fourier transforms of the transverse electric and mag¬ 
netic fields (taking into account that they depend only 
on the cca coordinate) 


Ek,i [t) 


Bk,i (t) 


1 


it), 




( 22 ) 



where i = 1,2, and further decompose the Fourier- 
transformed fields into the helical components Ek^R/L (t) 
and Bk^R/L with right- and left-handed helicities: 

BkM (t) = 2 {Bk,i (t) + B-k,i (t)) + 

+ ^ {Bk,2 (t) — S-fc 2 (t)) , 

2 i 

Bk,L {t) = — {Bk i (t) — B_k i (t)) + 

2 i 

+ 2 (^) + B-k,2 (t)) ■ (23) 

For electric fields, the definition of helical components 
is exactly the same. Again, here the term “helicity” 
refers to the direction of rotation of transverse electric 
and magnetic fields along the spatial direction of wave 
propagation (the X 3 axis in our setup). 

After such a decomposition, we calculate the energies 
of left- and right-handed helical electric and magnetic 
fields with a given wave number k as 

Ik.R/L (^) = \Bk,R/L it)f /2 -I- \B_k,R/L (0^/2, 

IkM/Lit) = \Ek,R/L{t)\^ /2+\E_k,R/L{t)\^ /2, (24) 

where k = and m = 0 ... [A 3 / 2 J now spans only the 
half of the discrete lattice momenta. Since the initial con¬ 
figuration of electromagnetic fields contains plane waves 
with (random) linear polarizations and equal energies, at 
t — 0 the energies Ik,R/L (t) of all left- and right-handed 
electromagnetic modes with 0 < k < are equal. 

We have found that for all simulations the energies 
^k R/L (^) exhibit quite large short-scale fluctuations with 
period of order At ~ 10 ... 100, which is smaller for short- 
wavelength modes and larger for long-wavelength ones. 
For illustration, see Fig. 7, where we plot the time de¬ 
pendence of Ik R (t) within a short initial period of time 
for simulation with parameter set No. 1. These oscil¬ 
lations indicate that the helical magnetic and electric 
fields represented by the basis (23) are not the eigen¬ 
states of the evolution process, which is in sharp contrast 
to the solution (4) of the anomalous Maxwell equations 
(2). We have explicitly checked that if the backreaction 
of fermions on the electromagnetic fields is neglected. 


FIG. 7: Time dependence of the energies (t) of right- 

handed components of electric field on a short time interval at 
the beginning of evolution for parameter set No. 1 {La = 200, 
n = 10, ra = 0.75, / = 0.2). The wave numbers are coded in 
color, from pure blue for the smallest nonzero value k — ^ 
(largest wavelength) to pure red for k = . 

these oscillations disappear and the energies (t) 

are constant in time for all values of k and for all polar¬ 
izations. This observation suggests that the short-scale 
oscillations might originate from the nontrivial depen¬ 
dence of fermionic current on the frequency, wave number 
and amplitude of electromagnetic field, which turns the 
solutions of the anomalous Maxwell equations (2) into 
waves with generic elliptic polarizations. 

Despite the short-scale fluctuations, we still find it use¬ 
ful to decompose our fields in the basis (23), since the 
corresponding electromagnetic modes carry definite he¬ 
licity and thus the energies of helical modes can be used 
to define, at least approximately, the helicity on the lat¬ 
tice. This definition is advantageous since direct lattice 
discretizations of the continuum formula "H ~ / dPxA ■ B 
are in general flawed by lattice artifacts. In order to ab¬ 
stract ourselves from the short-scale fluctuations, we de¬ 
fine the energies Ik,R/L (t) which are averaged over some 
finite time interval T: 

t+T/2 

^k’R/L — 'if j d,t {t ) . (25) 

t-T/2 

We have used the value T = 25, which is sufficient to 
remove practically all short-scale oscillations. 

On Figures 8 and 9 we separately illustrate the time 
dependence of the energies of the left-handed (on the left) 
and right-handed (on the right) helical magnetic and elec¬ 
tric fields with wave numbers k < for several most 
characteristic sets of simulation parameters. The wave 
numbers are coded in color, from pure blue for the small¬ 
est nonzero value k = (largest wavelength) to pure red 
for k = Semi-transparent colored regions show the 
range of short-scale oscillations of (t), and thick 

solid lines show the time dependence of the time-smeared 
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FIG. 8: Time dependence of the energies of helical magnetic fields for several selected sets of simnlation parameters. The 
wave numbers are coded in color, from pure blue for the smallest nonzero value k = (largest wavelength) to pure red for 
k = Semi-transparent colored regions show the range of short-scale oscillations of (t), and thick solid lines show 

the time dependence of the time-smeared energies (t) defined in (25). In black-and-white version pure blue and pure red 

correspond to black and light grey, respectively. Horizontal green (light grey) lines show the initial energies which are equal for 
all modes. Left-handed and right-handed modes are in the left and in the right columns, respectively. 
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FIG. 9: Time dependence of the energies of helical electric fields for several selected sets of simulation parameters. The 
wave numbers are coded in color, from pure blue for the smallest nonzero value k = (largest wavelength) to pure red for 
k = Semi-transparent colored regions show the range of short-scale oscillations of n/]^ (t), and thick solid lines show 
the time dependence of the time-smeared energies (t) defined in (25). In black-and-white version pure blue and pure red 

correspond to black and light grey, respectively. Horizontal green (light grey) lines show the initial energies which are equal for 
all modes. Left-handed and right-handed modes are in the left and in the right columns, respectively. 
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energies Ik,R/L {i) defined in (25). Horizontal green lines 
show the initial energies which are equal for all helical 
components of electric and magnetic fields. 

From Fig. 8 we see that in some simulations (param¬ 
eter sets No. 3-7) the energies of the helical compo¬ 
nents of magnetic field exhibit the expected signatures 
of the inverse cascade due to chiral plasma instability 
[3, 6 , 13, 46]. Namely, the energy of a single longest- 
wavelength right-handed helical mode rapidly grows at 
early times and reaches some saturation limit at late 
evolution time, whereas the energies of all other modes 
decrease with time. As expected from the anomalous 
Maxwell equations (2) with acME = I^a/ (27I'^), increas¬ 
ing HA by a factor of two (to ha = 1-5, parameter set 
No. 3, second row in Fig. 8 ) results in the growth of 
two right-handed modes. Comparing the data on Fig. 8 
and Fig. 5, we conclude that the growth of helical mag¬ 
netic fields is not necessarily accompanied by the decay 
of the axial charge, and vice versa (see also Table I for 
a summary of all simulations). Yet another observation 
which supports this conclusion is that in simulations on 
the smallest lattice, for which the axial charge exhibits 
most rapid decay, we have not found any signatures of 
the growing electromagnetic fields. Interestingly, increas¬ 
ing the value of ha and/or the number of initially ex¬ 
cited electromagnetic held modes also does not necessar¬ 
ily speed up the inverse cascade and the decay of Qa- 

An even more interesting picture emerges if we also 
consider the energies of the helical components of elec¬ 
tric Helds, shown on Fig. 9. It turns out that for some 
simulation parameters the long-wavelength helical com¬ 
ponents of the electric Held, rather than the magnetic 
Held, are enhanced during the evolution (parameter sets 
No. 1, 2, 8 ). For parameter set No. 3, both magnetic and 
electric Helds grow in time. It is remarkable that pre¬ 
cisely for these parameter sets the axial charge exhibits 
most rapid decay. It seems that both the growth of the 
helical electric Helds and the decay of the axial charge 
are triggered by sufficiently large initial amplitudes of 
electromagnetic Helds. Thus it seems that the roles of 
electric and magnetic Helds in the chiral plasma insta¬ 
bility scenario are essentially diHerent, in contrast to the 
simple solution (4) of the anomalous Maxwell equations. 
It is also interesting to note that we observe the maximal 
growth of long-wavelength helical electric Helds in simu¬ 
lations with a smaller value of Fermi velocity vp = 0.75 
(parameter set No. 8 ). Such a strong dependence on the 
Fermi velocity calls for a proper theoretical analysis. 

In order to quantify the net transfer of energy due to 
the inverse cascade, we follow [46] and introduce the mag¬ 
netic and electric correlation lengths (t) and (t) as 

it) 

^ ( 26 ) 

k 

where (t) = {t) + (t). The time depen¬ 

dence of (t) and {t), shown on Fig. 10, quantiHes 


the direction of the transfer of energy between short- and 
long-wavelength modes, ^b {t) and ^e (t) can be also 
thought of as the average wavelengths of magnetic and 
electric Helds at a given moment of time. The data shown 
on Fig. 10 indicates that ^e and ^b on average increase 
with time practically for all our simulations, thus provid¬ 
ing a more quantitative evidence for the inverse cascade. 
The growth is somewhat more pronounced for the elec¬ 
tric correlation length ^e, especially in simulations with 
larger initial amplitude / = 0 . 2 . At late evolution times, 
^E saturates at its upper bound ^e = A 3 equal to the lat¬ 
tice size. In contrast, the magnetic correlation length ^b 
exhibits rapid growth only at early times, and later seems 
to saturate at values smaller than L 3 . On Fig. 10 we do 
not show the data for the parameter set No. 9, since in 
this case we have found that only a single initially ex¬ 
cited mode strongly dominates the spectrum throughout 
the whole evolution process, and the quantities ^e and 
^B are trivially equal to L 3 up to some very small cor¬ 
rections. 

The effect of saturation of ^e,b at late evolution times 
prevents us from checking the universal late-time behav¬ 
ior ^E,B ~ Vi which follows from the scaling analysis of 
the anomalous Maxwell equations [13, 46], similarly to 
(21). It seems, however, that the late-time behavior of 
is more similar to Vi than that of ^b- 

So far almost all theoretical studies of the chiral plasma 
instability assume that the axial charge is distributed 
homogeneously in space and can be described by a 
coordinate-independent chiral chemical potential ha at 
all evolution times. The extension of the anomalous 
Maxwell equations (2) which allows to consider spatially 
inhomogeneous distributions of axial charge density has 
been constructed only recently in [27]. It is thus inter¬ 
esting to check how well the assumption of spatial ho¬ 
mogeneity of the axial charge density qax holds in our 
simulations. In order to quantify the spatial inhomo¬ 
geneity of qAx, we consider the space-averaged squared 
deviation of qax from its space-averaged value QaIV: 


O' [qa] = 



QA/vf. 


(27) 


The time dependence of the ratio (j[qA\ / {Qa/V) is 
shown on Fig. 11 for those sets of simulation parame¬ 
ters which exhibit axial charge decay (in particular, this 
Hxes / = 0.2). Since the Hamiltonian which we use to 
deHne the initial state of our simulations involves the 
initial spatially inhomogeneous conHguration Aq of vec¬ 
tor potential, even at the start of the evolution the axial 
charge density is slightly inhomogeneous, with deviations 
from mean value being of order of 5 %. As one can see 
from Fig. 11, at late evolution times the inhomogeneity 
of qxA tends to slightly increase, however, this increase 
is not dramatic and does not exceed 20%. This suggests 
that the approximation of spatially homogeneous axial 
charge distribution is not unreasonable even when the 
long-wavelength modes are strongly enhanced and domi¬ 
nate the evolution. For simulations with / = 0.05 which 
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FIG. 10: Magnetic and electric correlation lengths and (defined in (26)) as functions of time for different parameter sets. 
A smearing procedure similar to (25) was applied in order to suppress minor short-scale fluctuations in the data. 



FIG. 11: Time dependence of the standard deviation of the 
axial charge density a [qa] from its volume-averaged value 
Qa/V in simulations with / = 0.2. 


do not exhibit the decay of the axial charge, the inho¬ 
mogeneities of the axial charge density remain approxi¬ 
mately constant or even tend to decrease. 

An interesting question is also the net transfer of en¬ 
ergy between fermions and electromagnetic fields. As dis¬ 
cussed in Section II, in classical statistical field theory al¬ 
gorithm the total energy of fermions and electromagnetic 
fields is conserved up to the work performed by the ex¬ 
ternal current (see Fig. 14 in Appendix A for a numerical 
demonstration). Since in the simulations of chiral plasma 
instability discussed in this Section the external currents 
are absent, the transfer of energy can be characterized 
by the time dependence of the energy of electromagnetic 
field alone, which is illustrated on Fig. 12. We see that in 
almost all simulations the energy of electromagnetic field 
decreases or stays constant. The only exception is the 
simulation with parameter set No. 7 (n = 10, / = 0.05, 
fj-A = 1.5), for which the energy of electromagnetic field 


quickly increases by almost a factor of three at t > 6-10^. 
Analysis of power spectra suggests that this increase can 
be at least partly attributed to the enhancement of helical 
long-wavelength electric helds (similarly to the one ob¬ 
served for parameter set No. 3, see second row on Fig. 9). 
The decrease of electromagnetic field energy in all other 
simulations indicates that it might be not completely cor¬ 
rect to think of chiral instability as of a “discharge” of 
on excited state of Dirac sea into electromagnetic waves. 


V. CONCLUSIONS 

In this work, we have studied the real-time quan¬ 
tum evolution of chirally imbalanced Wilson-Dirac lattice 
fermions coupled to dynamical classical electromagnetic 
field within the classical statistical held theory approach. 
The quantum evolution of fermions was simulated ex¬ 
actly (up to small fully controlled errors originating from 
discretization of time). Our simulations of the chiral¬ 
ity pumping process, described by the volume-integrated 
anomaly equation (5), suggest that the effect of explicit 
chiral symmetry breaking due to the Wilson term in the 
lattice Dirac Hamiltonian is not very large. We hope 
therefore that our results can be confronted at least at 
the qualitative level with the theoretical predictions for 
continuum chiral fermions. 

We have considered both the generation of chirality im¬ 
balance in parallel electric and magnetic fields and the de¬ 
cay of initially present chirality imbalance at the expense 
of generating electromagnetic fields with nonzero helicity. 
We have observed that in general the backreaction of dy¬ 
namical electromagnetic fields prevents fermions from ac¬ 
quiring large chirality imbalance - either by suppression 
of the chirality pumping or by accelerating the decay of 
initially present chirality imbalance. The suppression of 
the chirality pumping process can be understood as the 
dynamical screening of the external electric field, sim¬ 
ilarly to what happens in the Schwinger pair creation 
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FIG. 12: Time dependence of the total energy of electromagnetic field in simulations with / = 0.2 (on the left) and / = 0.05 
(on the right). A smearing procedure similar to (25) was applied in order to suppress minor short-scale fluctuations in the data. 


process [36, 38]. 

In simulations with nonzero initial axial charge Qa we 
have also found a numerical evidence of the inverse cas¬ 
cade phenomenon due to the chiral plasma instability - 
that is, rapid growth of long-wavelength magnetic fields 
of definite helicity at early evolution times and the decay 
of all other magnetic field components. In some cases, he¬ 
lical electric fields were found to grow, even when mag¬ 
netic fields did not exhibit any enhancement. A sum¬ 
mary of our simulations given in Table I suggests that 
the growth (or at least the absence of decay) of long- 
wavelength helical electric fields is a necessary condition 
for the dynamical decay of the axial charge. The fact 
that the enhancement of helical electric fields is switched 
on only for sufficiently large initial amplitude of electro¬ 
magnetic field indicates that nonlinear responses such as 
the dynamical refringence [53] might be important for 
the evolution of chirally imbalanced plasma. 

We have observed the mechanism which eventually 
stops the growth of long-wavelength modes in our sim¬ 
ulations is not directly related to the decay of the ax¬ 
ial charge. This observation, together with quite differ¬ 
ent roles of electric and magnetic fields in the evolution 
process, suggests that the nontrivial momentum and fre¬ 
quency dependence of both the electric conductivity and 
the chiral magnetic conductivity might be important for 
the quantitative description of chiral plasma instability. 
On the other hand, our simulations also indicate that the 
approximation of spatially homogeneous axial charge dis¬ 
tribution, assumed in most theoretical considerations of 
anomalous Maxwell equations, is reasonably good even 
at late evolution times, when the instability has fully de¬ 
veloped and the growth of long-wavelength helical elec¬ 
tromagnetic fields has saturated. 

An interesting further development of our work would 
be to use chiral lattice fermions in the CSFT algorithm, 
with the possible choice of overlap Hamiltonian [54]. In 
this case, axial charge is conserved in the absence of elec¬ 
tromagnetic fields, and the effects of explicit chiral sym¬ 


metry breaking at high momenta should be absent. Such 
setup should be more relevant in the context of high- 
energy physics, where chiral symmetry tends to be exact 
at sufficiently high energies (at least at the level of the 
bare Lagrangian). Yet another interesting open question 
is the effect of the quantum fluctuations of the electro¬ 
magnetic field, which are encoded in the nontrivial initial 
density matrix. 
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Appendix A: Classical statistical field theory 
algorithm 

The starting point of our derivation of the CSFT al¬ 
gorithm is the general expression for the time-dependent 
expectation value of some quantum operator O: 

(O (t)) = Tr (poC (to, t) OtJ^ (to, t)) , (Al) 

where pq is the initial density matrix and the evolution 
operator tj (to,t) is the time-ordered exponent 

[>(to,t) =rexp (it'iF(t')j , (A2) 
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where the Planck constant is set to one by an appropri¬ 
ate choice of units. The Hamiltonian operator is defined 
by equations (6), (7) and (9) in Section If. We have al¬ 
lowed for an explicit time dependence of the Hamiltonian, 
for example, due to the time dependence of the external 
current Jx^i (t). The evolution operator (A2) can be ex¬ 
panded into a product of elementary evolution operators 
for small time step S = where should be much 
larger than any relevant energy scale in the system: 

U {t,to) = 

= lim (^e-^^Hito)s^-iH{to+s)s ^ ^ _ (^3) 

Let us now insert the decompositions of the identity op¬ 
erator ^ = n/ dAx^i \ Ax^i){Ax^i\ in the Hilbert space of 

x,i 

electromagnetic field between the infinitesimal factors as 
well as at the beginning and at the end of the product in 
(A3) in order to arrive at the path integral representation 
of the evolution operator (A2). We do this both for the 
forward and the backward evolution operators U {to,t) 
and in (Al). It is convenient to enumerate 

the gauge fields which enter identity decompositions in 
the forward evolution operators with the discrete lattice 
time variable t = 0 ... N, and in the backward branch - 
with T = 77 -I- 1... 2N + 1. The variable r is a discrete 
parametrization of the Keldysh contour going from to to t 
and back (see Fig. 13 for an illustration). Now we have to 
express the matrix elements {A'^\ \ of the 


elementary evolution operators in terms of the fields A'^ 
and In the derivation of the CSFT algorithm, it is 

most convenient to use the approximate expression 

(A^l |A^+i) « x 

xe (A4) 

for the forward evolution operators, and the different ap¬ 
proximate expression 


xe • (-^ 5 ) 

for the backward evolution operators. In the first expres¬ 
sion (A4), we order the electromagnetic field operators 
Ex^i and Ax^i in such a way that all the operators in 
the exponent containing Ax^i act on the vector {A'^\. In 
the second expression (A5), these operators act on the 
vector These approximations are both valid to 

order 0{5) and differ only in the terms of order 0(5^), 
hence being equivalent in the limit 5 —>■ 0. 

Using (A4) and (A5), we arrive at the path integral 
representation of the expectation value (O (t)), in which 
we integrate over the gauge fields living on the discretized 
Keldysh contour: 


{d{t)) = j dA° ... dA^^+^pEM [A°, A^^+Y X 

^ -^SHAA0] 

X ir p_p e ’ e ie ’’ 


X ... 


7-iV-l 


X . . . 


... X e 


- r iv+21 f T, (F^+A^+z 5 J:A^+^J 


(A6) 


It is important to stress that at this point we have used 
the path integral representation only for the bosonic 
fields, and the exponential factors e^'‘^FFlA ] (A6) 

are still operators in the fermionic many-body Hilbert 
space. Correspondingly, the trace in (A6) is taken 
over this Hilbert space. The operator of the observ¬ 
able 0 \_A^,A^'^Y = (A^|0|A'^+^) is also an oper¬ 
ator on the fermionic Hilbert space which depends on 
fields A^ and In deriving the above path inte¬ 

gral representation, we have made a simplifying assump- 



^ # ••• ^ ^ < 

t=2N+1 t=2N t=N+2 t=N+1 


FIG. 13: An illustration of the Schwinger-Keldysh contour 
with the discrete lattice time r. 
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tion [38] that the initial density matrix po factorizes into 
the direct product of the fermionic density matrix pp 
(which might in general be correlated with the initial 
state of the electromagnetic field) and the density matrix 
Pem of the electromagnetic field with matrix elements 
Pem T = {A^\pem\A^^'^^). While this as¬ 

sumption is certainly not valid, say, for the density ma¬ 
trix p = describing the thermal equilibrium state 

of the full Hamiltonian Hem + Hp, it is still justifiable in 
the case of almost classical dynamics of electromagnetic 
fields. 

At this point let us assume that the observable oper¬ 
ator O can be represented as a sum of the 

identity operator in fermionic Hilbert space (this sum¬ 
mand corresponds to purely bosonic observables) and of 
all possible fermionic bilinear operators: 

6 (A^, A^+i) = Ob (A^, A^+^) I -h 

x,y 

This form is sufficiently general to describe all the ob¬ 
servables which we consider in this work. Furthermore, 
let us assume that the fermionic density matrix pp can 
be represented as an exponent of some fermionic bilinear 
operator i7o = E '^v' 

pp = Z-i exp (-Hq/T^ , (A8) 

where T is some (perhaps fictitious) temperature. Note 
that in the case of evolution which starts from non¬ 
equilibrium state the operator Hq can be different from 
the Wilson-Dirac Hamiltonian Hp which governs the 
quantum evolution. For instance, the excited state with 
initial chiral imbalance considered in Section IV corre¬ 
sponds to the following form of Hq in the limit T —>■ 0: 

h^=h [A°] + PA^ |V'a)sign ((V'ahs \^a)) (V'al ,(A9) 

a 

where \ipa) are the eigenstates of the Wilson-Dirac 
Hamiltonian (7) with the initial gauge field Aq. It is 


obvious that for exactly chiral Dirac Hamiltonian with 
(V’al 75 this definition reduces to the form 

ho = h + pAlb- 

Now we are in the position to further simplify the trace 
over the many-body fermionic Hilbert space in (A6). To 
this end we use the identities 

Tr = det {l + , (AlO) 

Tr ... e^^Op^ = det (l -|- ... e^") x 

xTr (^{1 + e-^- Op^ , (All) 

where the operators Bi = ijjy and Op = 

x,y 

E i’lOx,yy are the fermionic bilinear operators, and 

x:,V 

the corresponding symbols without hats denote opera¬ 
tors on the single-particle fermionic Hilbert space with 
matrix elements and Opx,y Correspondingly, on 

the left hand side the traces are over the many-body 
fermionic Hilbert space, and the determinants and traces 
on the right hand side are on the single-particle fermionic 
Hilbert space. 

As yet another preliminary step in the derivation of 
the CSFT algorithm, let us also decompose the gauge 
fields on the forward and the backward branches of the 
Keldysh contour into the “classical” gauge field A” ^ and 
the “quantum” gauge field A” ^ as 

4” =4” -4- - P 

■^x,i ■^x,i ' 

= Ay-\Ay, t = 0...N. (A12) 

Relying on the assumptions (A7) and (A8) and us¬ 
ing the identities (AlO) and (All), one can rewrite the 
expression (A6) in terms of the operators on the single¬ 
particle fermionic Hilbert space and the variables A” ^ 
and A”_j: 


[d (t)) = Z ^ j dA° ... dA^ J dA° ...dA^ Pem (aP + ^.AP— ^ j exp ^Tr \tl{i u-e x 


X exp 


. N-l 
I 

d 


N-\ 


Y. Y [Kf - K!) (&- ■«;.) -i^YY +Y 7',.,, - 


r=0 x,i 


T—0 x,i 


X On A^ 


+ ^ + li' ( (l + O. - 4 " - ^ ) ) ) . (A 13 ) 


A^ 


A^ 


A^ 


where the field strength tensor F” is constructed from the “classical” component of the gauge field A” ^ exactly 
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in the same way as in (10) and we have introduced the 
unitary single-particle forward and backward evolution 
operators 




(A14) 


The path integral representation (A6) is exact in the 
limit A^ —>■ oo, (5 —?► 0 with fixed t = NS (up to the 
simplifying assumptions on the form of the observable 
operator O and the initial density matrix p), but is 
not suitable for numerical analysis. The key step in 
the derivation of the CSFT algorithm is to expand the 
fermion-induced effective action of electromagnetic field 
Sp = Tr In (l -I- in the first line of (A13) to 

the linear order in the “quantum” electromagnetic field 
A 


N 


Sf « 'S'j’U-.=o + Z^Z^ 


r=0 x,i 


dAl 


-Sf 


.(A15) 


.=0 


In order to calculate the first derivative oi Sf over A” j, 
we use the identities 


d 


dAl 


u+ 


Al , =0 


iS 

= - —M (0, r) j [A”] u (r, A^) 


dAl 


id 


.=0 


= (r, N) j [A”] (0, r) ,(A16) 


where we have introduced the single-particle operator of 
the conserved electric current 


jx,i [A] 


dh [A] 

dAx.i 


(A17) 


as well as the single-particle evolution operator in the 
background of the “classical” electromagnetic field A” 


u{ti,T 2 ) = T 2 > ri.(Al8) 

The identities (A16) are exact up to the order O (<5^), 
since in the derivatives of the forward and backward evo¬ 
lution operators we have used different orderings of the 
elementary evolution operator ] and the current 

operator j [A”]. 

Using (A16) and the relations u+(0, A^)|^^g = 

u{0,N) and U- (0, A^)|^^g = u'^ {0,N), after some sim¬ 
ple algebraic manipulations we can rewrite the derivative 
over A” ^ in (A15) as 


d 


dA. 


^{3h) = 


= -iSTr 


1 


1 -|- e^olT 


Al ,,=0 

u{0,T)jx,i [A”] u'^ (0,t) ) (A19) 


Now that our action (A15) is assumed to be lin¬ 
ear in the “quantum” field A” ^ for r = 1... A^ — 1, the 

quantum field A” j can be integrated out in a straight¬ 
forward way in the case of purely bosonic observables 
with Of = 0. The case of fermionic observables 
with nontrivial Of (a^ + A^ - is more sub¬ 
tle, since in the path integral representation (A13) the 
fermionic observable itself depends on A” ^ (via the fac¬ 
tor (l + ^ under the fermionic trace in 

the last line of (A13)). It is a common assumption in 
the derivation of the CSFT algorithm to neglect the 
A” j dependence of the fermionic observables (see e.g. 
[38]), which can be justified, e.g., if the relevant physi¬ 
cal processes involve large number of virtual fermionic 
particles. In this case one can argue that the expo¬ 
nent of the effective action Sf has a much stronger de¬ 
pendence on A” j than the observable. A heuristic ar¬ 
gument in favor of such assumption is that if one ne¬ 
glects the A” j dependence of the fermionic observables, 
the expectation values of all fermionic bilinear opera¬ 
tors take exactly the same form as the expectation val¬ 
ues of the electric current (A19) and the fermionic en¬ 
ergy (see equation (15) below). Since these quantities 
are related to the observables characterizing the classi¬ 
cal electromagnetic field via the inhomogeneous Maxwell 
equations and the energy conservation law, they are 
certainly also physical observables. While the A” de¬ 
pendence of the observable operator might still encode 
some interesting effects of the backreaction of measure¬ 
ments on the quantum evolution, taking it into ac¬ 
count would presumably lead to a significant compli¬ 
cation of the CSFT algorithm. For all these reasons, 
we also assume that the factor (l -|- in 

(A13) depends negligibly weakly on A” and replace it by 
(l -f M-i (0, N) (0, 

In order to integrate out the fields A° ^ and A^j at 
the endpoints of the Keldysh contour, it is convenient 
to introduce the Wigner transforms pem ^x,i) 

Oi ,2 (A^j, E^i) of the initial density matrix pem and the 
operators Of,b in (A7): 


Psm(a° + ^,A°-^) = 


= J dEx^^PEM (A°,J, Ex^i) 


(A20) 


(A21) 


where A” ^ is the “classical” electric field. We also note 
that the first sum over r in the second line of (A13) can 
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be rewritten as 

T—O x,i 
N-1 

= “ E E ^ ^ ~ 2-4^,i) + 

T — 1 x,i 


+J 2 (^5 - <r') - - ^°..)) (A22) 

x,i 

Finally, we are ready to integrate out the “quantum” 
electromagnetic field A^. which leads to the following 
expression for the expectation value {O (t)): 


[6 {t)) = j dE^ dE^ J dA°... dA^pEM {A°, E°) x 


xS 


E° -;- STZ° 


N-1 


^r+l _ 2 ^t 


(57^" 


E^ - 


A^ - 


X +Tr + \ {0, N) Oi {A^, E^) u^O, , (A23) 


where 

= Jh + (jh) + E ’ (A24) 

j 

and {jli) is the expectation value of the electric cur¬ 
rent defined as in (A19). We note that the normalization 
factor Z~^ in (A8) and (A13) is cancelled by the zeroth- 
order term of the expansion (A15). 

From the explicit expression (A19) for the fermionic 
electric current (^) one can immediately see that it 

depends only on the classical electromagnetic field A” ^ 
with r' < T. Therefore the delta-functions in the in¬ 
tegral (A23) can be regarded as the constraints on the 
deterministic evolution of the classical electromagnetic 
field j interacting with the quantum fermionic field. 
To make this more obvious, we can rewrite the chain of 
(5-functions in (A23) as 


all the other delta functions. Integrating out A^, we 
remove the second delta-functions and replace A^ by 
[A°,£'°] = (A^ [A°,A°] ,A°,(j^)). We can re¬ 

peat this process for all r up to A^ — 1, each time ex¬ 
pressing A” in terms of the initial values A° and E^ and 
the functionals A^ with r' < r. It is important that in 
such a sequential integration, the integrand A"^ always 
enters the argument of the delta-function being removed 
linearly. Therefore despite the nonlinearity of the chain 
of evolution equations, such intermediate integrations do 
not produce any nontrivial Jacobian. To our knowledge, 
the absence of Jacobian in the integration measure in the 
CSFT algorithm so far has only been demonstrated for 
scalar field theory [-55, 56]. It is nice to see here explicitly 
its absence for lattice gauge theory coupled to fermions. 

After integrating out all the intermediate fields A"^ with 
T = 1... A^ — 1, we are left with the following form of 
equation (A23) 


J [A^ - A^ [a°,f;°]] X 

N 

xJ]<5[A"-A" [A"-\A"-2,(r-')]] X 

t=2 


xS 



6 


A^ [A°, A°] =A° + S{E°- (57^°), 
A^ [A”-i,A”-2,(r-i)] = 
= 2A^-i - A”-2 - (527^^-^ 


(A25) 


{6 {t))= J dA° dE° J dA^ dE^pEM (A°, E°) x 
x6 [A^ - A^ [A°, A°]] J [E^ [A°, A°]] X 

X (Oo {A^,E^) + 

+T" { l + lho/T ^ (O’(0, Af))) a26) 


where 


[A°,E°] 


an ^a°,E°] - a^-1 [A°,E°] 
~5 


(A27) 


where the last definition is for t = 2 ... N. From this 
expression one can see that one can sequentially inte¬ 
grate out the fields A'’’ with r = 1... A^ — 1 and express 
the fields A-^, E^ in terms of the initial values A°, 
E^. Namely, integrating out the field A^ first, we re¬ 
move the first delta-function in the product in (A25) 
and replace A^ by A^ [A°,A°] in the arguments of 


and we have expressed the functionals A'^ and A'^”^ in 
terms of the initial values A°, E^ of the vector gauge 
potential and the electric field. In this expression, it 
is straightforward to integrate out A^ and E^, which 
amounts to substituting A'^ [A°,A°] and [A^,!?^] 
in place of A^ and E^ in the Wigner transforms of the 
observable operators Ob and Of- 


















20 


To summarize, the CSFT algorithm amounts to a si¬ 
multaneous time evolution of the classical electromag¬ 
netic helds, described by the vector potential A^. ^ and 
the electric field and the quantum fermionic fields, 
described by single-particle evolution operator u (0, r) in 
(A18). The discrete equations which govern this evolu¬ 
tions (arguments of the delta functions in (A23)) have a 
well-defined continuum limit 5 —>■ 0, N —>■ oo with fixed 
t = NS, which is given by the equations (11), (12) and 
(13) in the main text. To simplify the notation, in the 
main part of the text we omit the bar over the classi¬ 
cal components of the gauge field and denote them as 
Ax,i it) = Ex^i (t) = Fx^ij (t) = 

In practice, however, the numerical solution of the con¬ 
tinuum equations (11), (12) and (13) should necessar¬ 
ily involve some discretization of time. While the sim¬ 
ple discretization of the Keldysh contour used in the 
above derivation can be in principle used for numerical 
solution at sufhciently small 5, for a given finite 5 one 
can construct different, more advanced discretizations 
which would reduce discretization errors, thus improving 
the conservation of energy (15) and making the single¬ 
particle evolution operator u (0, r) numerically closer to 
a unitary matrix. 

In this work we follow [35, 38] and use the leapfrog 
evolution scheme for the single-particle evolution opera¬ 
tor vS' = u(0,r), which significantly improves the con¬ 
servation of the unitarity condition u (0, t) (0, r) = 1 
at finite discrete time step 5: 

- iSh [A^] u^, t = 

v} = — iSh [A°] = 1. (A28) 

In practice it is convenient to work with the components 
of u'^ in the basis of eigenstates of the initial single¬ 
particle Hamiltonian h [A°]. In particular, if transla¬ 
tional invariance along some of the space directions is pre¬ 
served during the evolution, remains block-diagonal 
in the basis of plane waves propagating along these di¬ 
rections. This block-diagonal structure can be used to 
greatly reduce the number of independent components 
of which enter the equations (A28). We have used 
translational invariance in two out of three spatial lat¬ 
tice directions to speed up the evolution algorithm on 
large lattices with sizes up to 200 x 40 x 40, assuming 
translational invariance in two out of three spatial direc¬ 
tions. 

For the evolution of electromagnetic field we use the 
equations which directly follow from (A23): 


tt’T+I rpr 

- ^xA 




4'^+! _ AT 

'^x,i -^x,! _ E^r+1 

~ ^ ^xA ’ 


■ — AP 


_ 77 U 

-^xA 



t/10® 


FIG. 14: Time dependence of the energies of fermions and 
electromagnetic fields and their total for simulations on 200 x 
20 X 20 lattice with n = 10, / = 0.2, fiA = 0.75, wf = 1 
(parameter set No. 1 in Table I). 



In our simulations, we use the value <5 = 0.05. We have 
checked that decreasing <5 down to 0.02 does not change 
our results up to some small unimportant fluctuations. 

In principle, leapfrog-type time discretization (A28) al¬ 
lows the existence of fermionic doublers in time direc¬ 
tion - that is, the symmetric finite differences in (A28) 
are zero if the mode functions oscillate as (—1)^. Such 
doubler modes correspond to another flavour of Dirac 
fermions with an opposite signature of the 75 matrix. 
Thus if such modes are excited, they can also contribute 
to the anomaly equation (5) and effectively decrease the 
anomaly coefficient, or lead to the decay of the initial 
value of the axial charge [33] . In order to check whether 
fermionic modes with such high frequencies are excited 
we have calculated the average norm of forward finite 
differences of as j^Tr 

Since the size of the single-particle Hilbert space is equal 
to AV this quantity should be of order of <5^ if u'^ are 
smooth functions of r. On the other hand, doubler modes 
with u'^ ~ (~1)’' yield the contribution of order of unity. 
In our simulations we have checked that the above norm 
remains of order of 10“^ for all evolution times and does 
not exhibit any tendency to grow. This suggests that the 
doubler modes remain practically unexcited during the 
evolution. 

Another important characteristic of the discretization 
of the evolution equations (11) and (13) is the precision 
with which the conservation of energy (15) holds. For 
the leap-frog equations (A28) and (A29) the total energy 
of electromagnetic fields and fermions is conserved only 
approximately, up to the terms of order of S^. In order to 
illustrate the conservation of energy in our simulations, 
on Fig. 14 we show the time dependence of the fermionic 
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energy {Hf), the energy Hem of electromagnetic fields 
and their total. One can see that while both {Hf) and 
Hem change quite strongly during the evolution, their 
sum is conserved with a very good precision, which again 
shows that the time step <5 = 0.05 is small enough. 

The discrete evolution equations (A28) and (A29) are 
ideally suited for parallelization on multi-node comput¬ 
ers. Indeed, the largest amount of computer time is 
required to solve the evolution equation (A28) for the 
(4F) X {AV) matrix . Taking into account that only 
half of the single-particle fermionic states are filled at 
zero temperature, this size can be reduced by a factor 
of two down to {2V) x (dF). The evolution of elec¬ 
tromagnetic field (A29) requires only the total electric 
current summed over all fermionic modes, and is com¬ 
putationally very cheap. Thus it is natural to distribute 
the rows of the matrix over multiple nodes. On each 
node, one performs the elementary evolution step (A28) 
for the rows attributed to this node and calculates the 
partial traces of the electric current and other fermionic 
bilinear operators over these rows. The results are sent 
to the master node, which calculates the total current 
and performs the evolution of the electromagnetic field 


according to ( A29) . Since the amount of data transferred 
by network from each slave node to the master node is sig¬ 
nificantly smaller than the amount of data stored at each 
slave mode (for realistic lattice sizes and node numbers, 
the number of rows of u'^ per node is large), the speed of 
the algorithm scales practically linearly with the number 
of slave nodes. 

Such parallelization also solves the problem with very 
large amount of RAM memory required to store the ma¬ 
trix (~ 36 Gb for the 20 x 20 x 20 lattice when one 
uses 8-byte double accuracy numbers for all fields), which 
is simply split over different nodes. In order to further 
decrease the required RAM size, we use the d-byte float 
numbers to store . We have explicitly checked that the 
reduction from double to float real numbers practically 
does not affect our results. 

The extensive parallelization also allows us to avoid the 
stochastic summation over all modes [35-38], which in¬ 
troduces additional statistical noise in the results and can 
therefore significantly affect potentially unstable evolu¬ 
tion which we study in this work. Instead, we perform ex¬ 
act summation over all initially occupied fermionic states. 
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